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NUMERICAL  SIMULATION  OF  THE  COMPRESSIBLE 
ORSZAG-TANG  VORTEX 

I.  INTRODUCTION 

The  theoretical  and  numerical  study  of  turbulence  has  concentrated  until  recently  on 
incompressible  media.  The  motivation  for  this  concentration  of  effort  is  the  comparative 
simplicity  of  the  incompressible  case — the  system  of  governing  equations  for  the  incom¬ 
pressible  case  is  of  lower  order  than  the  compressible  system,  and  an  important  constraint 
(solenoidality)  links  different  components  of  the  velocity.  In  addition,  incompressible  fluids 
have  ideal  invariants  which  facilitate  analysis  of  the  states  to  which  turbulent,  dissipative 
fluids  relax  (e.g.,  [1]).  Unfortunately  many  turbulent  media,  such  as  nonconducting  neu¬ 
tral  gases  and  space  plasmas,  axe  compressible,  and  the  relevance  of  the  incompressible 
results  to  their  behavior  is  questionable.  Marsch  and  Mangeney  [2],  for  example,  recognized 
this  uncertainty  regarding  the  treatment  of  turbulence  in  the  solar  wind  and  interplane¬ 
tary  medium,  and  derived  a  compressible  formulation  of  magnetohydrodynamics  (MHD) 
in  terms  of  Elsasser  variables  for  application  to  such  problems.  Our  purpose  here  is  to 
begin  the  systematic  development  and  study  of  numerical  data  to  determine  the  effects 
of  compressibility  on  MHD  turbulence.  The  data  also  provides  a  basis  for  interpreting 
theoretical  treatments  and  devising  a  more  complete  theory.  In  this  paper  we  examine  a 
set  of  initial  conditions  that  are  defined  analytically  for  both  incompressible  and  compress¬ 
ible  media.  By  comparing  the  results  of  numerical  simulations  of  the  compressible  case 
with  incompressible  simulations,  we  identify  some  of  the  effects  which  are  introduced  by 
compressibility. 

To  place  the  present  paper  in  perspective,  we  note  that  at  least  three  other  major 
categories  of  compressible  turbulent  phenomena  have  been  pursued.  The  most  obvious 
is  the  interaction  of  a  sound  or  shock  wave  with  a  turbulent  field  or  vortices  (see,  for 
example,  papers  by  Ting  [3],  Zang  et  a/.[4],  and  Ribner  [5]).  Here  the  emphasis  is  on 
the  amplification  of  turbulence  or  the  production  of  sound.  Closely  related  is  the  direct 
production  of  sound  by  a  turbulent  field,  which  was  studied  in  detail  by  Lighthill  [6  , 
7].  Third  is  the  production  of  vorticity  via  compressible  effects  in  an  initially  irrotational 
flow.  This  can  occur,  for  example,  through  the  baroclinic  source  term  in  the  equation  for 
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the  evolution  of  the  vorticity,  as  shown  by  Picone  et  al.  [8,9,  10].  The  vorticity  can 
subsequently  evolve  into  turbulence  by  interaction  of  resonant  harmonics  [11]. 

This  paper  treats  a  fourth  category,  the  evolution  of  an  unbounded  compressible  turbu¬ 
lent  field.  Here  the  medium  may  interact  with  itself  through  internally  generated  acoustic 
or  magnetoacoustic  waves.  The  property  of  compressibility  also  leads  to  alterations  in  the 
morphology  and  description  of  the  medium  exclusive  of  such  waves.  Thus  our  concern  is 
with  the  dynamics  of  “compressible”  turbulence  itself,  and  a  complete  treatment  requires 
a  fully  compressible  representation  of  the  entire  flow  field,  including  the  rotational  compo¬ 
nent.  To  achieve  a  useful  description  of  compressible  MHD  turbulence,  we  use  numerical 
experiments  on  analogous  incompressible  and  compressible  problems.  This  is  different  from 
the  first  two  research  areas  mentioned  above,  in  which  the  standard  approach  has  been  to 
treat  the  turbulent  portion  of  the  fluid  as  incompressible  and  separable  from  compressible 
wave  phenomena  such  as  shocks  and  sound  waves.  Such  a  separation  of  the  incompressible 
and  compressible  parts  of  the  flow  is  consistent  with  linearization  of  the  governing  equa¬ 
tions  (e.g.,  [12]).  Our  approach,  however,  requires  the  solution  of  the  full  set  of  nonlinear 
governing  partial  differential  equations  because  we  treat  the  entire  fluid  as  compressible. 

An  accurate  solution  of  the  nonlinear  equations  governing  a  compressible  turbulent 
field  includes  phenomena  associated  with  acoustic  or  shock  waves,  but,  in  our  case,  not 
as  external  sources  of  excitation  nor  as  actions  of  the  turbulent  field  on  external  regions. 
Thus,  the  turbulence  can  produce  acoustic  and  shock  waves  which  interact  with  each 
other  and  with  the  turbulent  field  itself.  In  addition  to  altering  the  turbulent  field,  these 
interactions  can  generate  vorticity  through  the  baroclinic  source  term  in  the  vorticity 
evolution  equation  [9].  Such  effects  are  thus  important  to  our  investigations  but  are  not 
the  whole  story,  as  confirmed  by  the  differences  between  incompressible  solutions  and 
homentropic,  polytropic  compressible  solutions  for  a  steady  state  vortex  in  an  unbounded 
medium  [13].  These  differences  exist  even  though  the  velocity  field  for  the  steady  state 
vortex  is  solenoidal  in  both  cases. 

This  paper  deals  with  initially  average  subsonic  flows  in  two-dimensional  compressible 
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magnetofluids  in  the  average  Mach  number  regime  from  0.2-0. 6  (The  limitation  to  Mach 
numbers  of  about  0.6  or  less  is  due  to  the  way  in  which  we  specify  the  fluctuating  pressure, 
as  described  below).  For  the  magnetofluid  there  are  also  compressible  alterations  to  the 
magnetic  source  terms  in  the  vorticity  equation,  which  indicate  the  possibility  of  vorticity 
production  by  the  interaction  of  magnetoacoustic  waves  and  MHD  shock  waves  with  the 
turbulent  field.  As  in  the  incompressible  case,  the  magnetic  field  is  still  constrained  to  be 
solenoidal.  Note,  however,  that  the  absence  of  the  velocity  field  solenoidality  constraint 
also  leads  to  the  presence  of  a  compression  term  in  the  magnetic  induction  equation.  This 
variation  in  mode  coupling  between  the  velocity  field  and  magnetic  field  implies  that  the 
magnetic  field  can  evolve  differently  in  the  compressible  medium;  for  example,  magnetic 
flux  can  accumulate  in  compression  regions. 

A  noteworthy  initial  approach  to  the  theory  of  compressible  MHD  turbulence  by  Mont¬ 
gomery  et  al.  [14]  expresses  the  fluctuating  density  as  a  function  of  the  pressure  derived 
through  the  usual  incompressible  formulation.  This  assumes  that  incompressible  flow  is  a 
well  defined  limit  of  compressible  flow  and  seeks  to  determine  the  spectral  index  of  density 
fluctuations  for  “nearly  incompressible”  MHD  turbulence.  Shebalin  and  Montgomery  have 
investigated  this  model  numerically  [15].  As  mentioned  above,  Marsch  and  Mangeney  [2] 
have  developed  a  theoretical  framework  which  extends  methods  used  in  the  incompress¬ 
ible  MHD  regime.  For  the  hydrodynamic  case,  Moiseev  et  al.[16]  have  used  similarity 
and  renormalization  group  arguments  to  estimate  the  deviations  of  subsonic  compressible 
turbulence  from  the  Kolmogorov  spectrum  in  a  noncharged  system  driven  by  a  Gaussian 
external  force. 

Recent  computational  studies  of  fully  compressible  turbulence  for  noncharged  fluids 
by  Passot  and  Pouquet  [17]  and  Erlebacher  et  <z7.  [18]  and  for  MHD  by  Dahlburg  et  al.  [19] 
have  used  random  initial  conditions  with  a  range  of  spatial  scales.  All  of  these  calculations 
used  spectral  methods  for  their  accuracy  and  cost-effectiveness  in  studying  turbulent  flows. 
Earlier  studies  of  compressible  MHD  turbulence  employed  finite  difference  algorithms  and 
initial  conditions  consisting  of  at  most  a  few  modes  [20,  21,  22]. 
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None  of  the  above  MHD  studies  directly  compared  compressible  and  incompressible 
simulations.  Erlebacher  et  al  (1987)  have  compared  compressible  and  incompressible 
time  histories  of  various  statistical  quantities  in  hydrodynamic  turbulence  [18].  Passot  and 
Pouquet  separated  the  solenoidal  and  nonsolenoidal  components  of  the  velocity  and  specific 
kinetic  energy  fields  in  order  to  accomplish  this  comparison  for  a  dissipative  fluid  obeying 
the  Navier-Stokes  equations  for  a  perfect  gas.  As  indicated  above,  however,  solenoidality  of 
the  velocity  field  is  not  a  sufficient  condition  for  incompressibility  of  the  medium,  so  that 
their  method  does  not  rigorously  separate  compressible  and  incompressible  phenomena. 
In  this  paper,  therefore,  we  take  a  somewhat  different  approach,  viz.,  we  compare  the 
evolution  of  compressible  and  incompressible  magnetofluids  starting  from  the  same  initial 
conditions.  We  perform  simulations  at  several  subsonic  initial  average  Mach  numbers  to 
determine  the  effects  of  increasing  compressibility.  We  use  the  nonrandom,  periodic  initial 
conditions  of  Orszag  and  Tang  [23].  The  initial  velocity  and  magnetic  fields  each  contain 
only  one  mode,  although  not  the  same  one.  The  level  of  compressibility  is  determined  by 
the  initial  mean  thermal  pressure  magnitude.  The  model  contains  an  incompressible  analog 
that  is  well  known  ( e.g [23],  [24],  [25],  [26],  [27]  [28]),  providing  an  excellent  baseline  case 
for  the  study  of  compressible  MHD  turbulence.  These  initial  conditions  thus  permit  a 
more  direct  (less  statistical)  quantitative  and  structural  comparison  of  turbulent  evolution 
within  compressible  and  incompressible  media  and  contains  most  of  the  significant  features 
of  MHD  turbulence,  including  dissipation  of  magnetic  and  kinetic  energy,  reconnection, 
formation  of  high  density  jets,  selective  decay,  dynamic  alignment,  and  the  emergence 
and  manifestations  of  small  scale  structure  (e.g.,  [29],  [30],  [31],  [32],  [33]).  A  preliminary 
version  of  this  work  has  appeared  elsewhere  [34], 

The  paper  begins  with  \  formulation  of  the  problem,  followed  by  a  section  describing 
our  numerical  methods.  Because  relatively  few  numerical  treatments  of  compressible  MHD 
turbulence  have  been  performed,  the  optimum  numerical  method  and  the  appropriate 
equations  to  integrate  are  issues  in  themselves.  In  the  course  of  our  discussion,  we  review 
formulations  and  algorithms  that  are  reasonable  candidates  for  consideration  and  have 
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been  employed  for  similar  purposes  or  physical  regimes.  Section  III  presents  our  numerical 
results,  including  the  significant  differences  between  the  compressible  and  incompressible 
analogs,  for  a  number  of  turbulent  MHD  phenomena.  Section  IV  presents  our  conclusions 
and  a  description  of  related  work  in  progress, 

II.  NUMERICAL  SIMULATIONS 
A.  Formulation 

We  start  with  the  nonlinear  partial  differential  equations  which  govern  the  behavior  of 
a  two-dimensional,  compressible,  dissipative  magnetofluid,  written  here  in  a  dimensionless 
form: 


dp 

dt 


Z  =  -V.(pv), 


d_{pv) 
dt 


— '  =  _  v  • 


pvv  -  BB  4-  \{p  +  |B|2)I  -  j-r] 
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supplemented  by  an  equation  of  state, 
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The  notation  in  Eqs.(l)  and  (2)  is  as  follows:  p  is  the  mass  density,  v  =  (vx,vv,0)  is 
the  flow  velocity,  p  is  the  mechanical  pressure,  I  is  the  unit  dyad,  B  =  (3z,By,0)  is  the 
magnetic  induction  field,  U  is  the  internal  energy  density,  E  is  the  total  energy  density, 
given  in  dimensionless  units  by 

E{x,  t)  =  p\x\2  +  |B|2  4-  U,  (3) 

T  =  p/ p  is  the  dimensionless  temperature,  r  =  -f  d,vj)  —  |V  •  v<5,y]  e,  e,  is  the 

viscous  stress  tensor,  and  7  =  5/3  is  the  ratio  of  specific  heats.  The  thermal  conductivity 
(«),  magnetic  resistivity  (7),  and  viscosity  (p)  are  constant  and  uniform,  and  we  assume 
that  the  Stokes  relation  holds,  so  that  the  bulk  viscosity  is  zero  [35].  The  important 
dimensionless  numbers  are  the  viscous  Lundquist  number  Sv  =  p0VaL0/p  [36],  the  resistive 
Lundquist  number  Sr  =  VaL0/t /,  the  Prandtl  number  Pr  =  cpp./n,  and  the  Alfven  number 
A  =  (YaI^ o2)*-  these  definitions,  po  is  a  characteristic  density,  VA  is  the  Alfven 

speed,  Lo  is  a  characteristic  length  equal  to  the  reciprocal  of  the  minimum  wavenumber 
kmtn  =  2n/L  ( L  =  system  length),  cp  is  the  specific  heat  at  constant  pressure,  and  V0  is  a 
characteristic  flow  speed.  The  thermodynamic  normalization  sets 

£0  =  Po  =  Bl/Sn  =  p0Vl/2  .  (4) 

The  characteristic  speed  of  sound  is  then  C|0  =  7p0/^o-  Time  ( t )  is  measured  in  units 
of  the  Alfven  transit  time  Lq/Va ■  For  the  runs  reported  here,  we  set  Sv  =  Sr  = 
50,  100,  or  200  and  Pr  =  A  =  1.  All  of  our  numerical  simulations  assume  unit  mag¬ 
netic  Prandtl  number,  i.e.,  p/(p0ri)  =  1,  so  we  subsequently  consider  only  the  parameter 
5  =  =  5r. 

The  reader  should  note  that  formulations  of  the  numerical  problem  other  than  Eqs. 
(1)  are  presently  being  used  to  study  compressible  MHD  phenomena.  Schnack  et  ai. [37]  use 
only  equations  (16),  (lc)  and  (le)  with  the  assumption  of  constant  mass  density  and  zero 
mechanical  pressure  to  model  the  Reversed- Field  Pinch  for  magnetic  fusion.  This  is  a  zero- 
/?  approximation  to  a  very  low  /?  plasma,  where  is  the  ratio  of  the  mechanical  pressure 
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to  the  magnetic  pressure.  Mikic  et  a/.[38]  also  have  used  ohis  approximation  to  model  the 
solar  corona  (but  note  that  they  replace  the  Eqs.  (lc)  and  (le)  with  a  magnetic  vector 
potential  equation).  In  addition,  they  assume  that  the  convection  term  in  the  equation  of 
motion  is  negligible.  The  applicability  of  a  zero-/?  model  in  these  two  cases  is  debatable. 
More  importantly,  mass  density  fluctuations  and  gradients  in  mechanical  pressure  can  play 
as  great  a  role  in  MHD  turbulence  as  they  do  in  nonmagnetized  fluids. 

More  common  in  treatments  of  compressible  flows  is  the  assumption  of  a  homentropic, 
polytropic  medium  [22],  [39]  which  satisfies  the  equation 

p  =  Q  p1  ,  Q  —  constant  .  (5) 

This  is  often  called  the  “polytropic  equation  of  state.”  Technically  speaking,  however,  a 
poly  tropic  fluid  need  not  have  constant  entropy,  or  equivalently,  Q  need  not  be  constant. 

Assuming  that  the  fluid  is  homentropic,  Eq.  (5)  replaces  Eq.  (Id)  for  conservation 
of  total  energy,  and  the  system  to  be  integrated  is  then  Eqs.(la)-(lc)  and  (le).  Note 
that  constant  entropy  implies  that  the  dissipative  terms  are  zero.  Unfortunately,  such  sys¬ 
tems  of  hyperbolic  equations  will  likely  develop  magnetoacoustic  waves  and  discontinuities 
(shocks)  in  a  finite  amount  of  time  when  the  flow  is  unsteady.  The  fact  that  the  entropy  is 
constant  and  uniform  throughout  the  fluid  means  that  the  shock  solutions  of  this  reduced 
system  of  equations  cannot  satisfy  the  Rankine-Hugoniot  jump  conditions,  except  for  very 
weak  shocks  [39].  Equation  (5)  also  precludes  the  generation  of  vorticity  through  the  baro- 
clinic  term  (^V/)  x  Vp)  in  the  vorticity  evolution  equation  [9].  Thus  the  assumption 
of  a  homentropic,  polytropic  medium  cannot  account  for  several  important  factors  in  the 
formation  and  evolution  of  compressible  MHD  turbulence.  Our  formulation  (Eqs.(l)  and 
(2))  provides  for  all  of  these  effects  and  for  strict  conservation  of  mass,  momentum,  and 
energy. 

The  use  of  Eq.(5)  instead  of  Eq.(ld)  also  has  important  consequences  for  numerical 
modeling.  In  an  ideal  fluid,  the  development  of  discontinuities  will  result  in  unstable  nu¬ 
merical  solutions  if  a  high  order  explicit  integration  scheme  is  used.  The  computer  model 
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will  therefore  "blow  up.”  The  addition  of  sufficiently  large  dissipative  terms,  either  phys¬ 
ical  or  numerical  ( e.g .,  a  viscous  term  in  the  equation  of  motion),  can  result  in  smooth 
mathematical  and  numerical  solutions  with  steep  finite  gradients  instead  of  shock-like  dis¬ 
continuities  f  'O].  When  Eq.  (5)  is  used,  such  terms  are  inconsistent  with  a  constant 
entropy,  as  mentioned  above.  In  addition,  if  this  “reduced”  system  of  equations  is  dissi¬ 
pative,  then  total  energy  is  not  conserved,  similar  to  the  case  of  incompressible  how.  The 
shocks  would  be  nonphysical  since  the  Rankine-Hugoniot  conditions  cannot  be  satisfied. 
The  size  of  the  coefficient  of  viscosity  that  is  necessary  to  guarantee  numerical  stability 
fox  an  explicit  algorithm  varies  with  the  situation.  A  fully  implicit  or  semi-implicit  inte¬ 
gration  scheme  can  ensure  that  finite  acoustic  or  magnetoacoustic  waves  are  filtered  out  of 
the  solution  [18,  41].  This  should  provide  a  stable  solution,  although  such  models  cannot 
represent  the  effects  of  finite  amplitude  waves  in  a  turbulent  compressible  fluid. 

B.  Initial  Conditions 

For  initial  conditions  we  use  the  vortex  system  of  Orszag  and  Tang  (1979): 

p(x,y,f  =  0)  =  1,  (6a) 

v(x,  y,  t  =  0)  =  —  sin  y  ez  -fi  sini  ey,  (6b) 

and 

B(x,  y,  t  =  0)  =  —  siny  er  +  sin  2x  ey,  (6c) 

whore  ex  and  ey  are  unit  vectors  in  the  x  and  y  directions  respectively.  The  initial  me¬ 
chanical  pressure,  p(t  =  0),  is  decomposed  into  a  mean  part,  (p(t  =  0)),  and  a  fluctuating 
part,  pj(t  —  0).  We  set  pj(t  =  0)  equal  to  the  appropriate  incompressible  pressure  distri¬ 
bution,  which  is  obtained  by  applying  a  divergence  operator  to  the  incompressible  MHD 
equation  of  motion  (Eq.  (lb')  below)  and  then  inverting  the  resulting  Poisson  equation. 
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This  specification  of  p/(t  =  0)  puts  a  physical  lower  limit  on  {p(t  —  0)),  i.e.,  we  must  select 
(p(t  =  0)}  such  that  p(x,y,t  =  0)  =  (p(x,y,t  =  0))  +  p/(x,y,t  =  0)  >  0  everywhere  in  the 
system.  Figure  1  shows  the  initial  form  of  the  flow  field  in  the  form  of  streamline  plots  of 
the  velocity  and  magnetic  field. 

Notice  that,  for  a  polytropic,  homentropic,  compressible  fluid  (Eq.(5)),  our  initial 
conditions  of  a  spatially  varying  mechanical  pressure  and  a  spatially  constant  mass  density 
would  contradict  Eq.(5).  However,  in  our  formulation,  the  temperature  varies  spatially 
with  the  mechanical  pressure,  so  that  the  mass  density  can  in  fact  be  constant.  Because  of 
this,  any  steady  state  solenoidal  velocity  field  that  can  exist  in  an  incompressible  medium 
has  a  direct  analog  in  a  compressible  medium.  An  alternative  to  our  initial  conditions 
would  be  to  use  Eq.  (5)  to  get  an  approximately  consistent  mass  density  field  rather  than 
Eq.  (la);  however,  the  mapping  onto  the  initial  conditions  of  the  incompressible  solution 
would  then  be  altered.  We  have  also  used  a  flat  initial  pressure  profile  along  with  Eqs. 
(6a)-(6c)  and  have  found  that  the  early  behavior  ( t  <  1.0)  is  somewhat  different  in  that, 
e.g.,  growth  of  correlation  is  more  rapid.  This  is  caused  by  the  small  expansion  of  the 
vortices  necessary  to  set  up  a  pressure  gradient  in  response  to  the  centrifugal  forces  that 
are  present.  In  the  calculations  of  Dahlburg  et  af.[19],  the  pressure  and  mass  density  were 
flat  and  broadband  noise  was  used  to  construct  the  magnetic  field.  The  initial  velocity 
field  resulted  from  rotating  the  magnetic  field  through  a  fixed  angle  consistent  with  the 
desired  value  of  the  initial  correlation  £.  Dahlburg  et  al.  found  that  again  the  evolution  of 
global  integrals  was  smooth  after  an  initial  rapid  jump  in  correlation. 

Returning  to  the  definition  of  initial  conditions,  the  normalization  in  Eq.  (4)  gives 
the  local  Mach  number  M.  as 


M2 


2pjvf 
7  P 


(7a) 


To  obtain  a  characteristic”  Mach  number  M  at  time  t  =  0  for  the  flows  in  our  calculations, 
we  substitute  p  =  (p(t  =  0)),  p  =  (p(t  =  0)),  and  |v|  =  \J (v2(f  =  0))  into  Eq.  (7a).  We 
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then  choose  initial  values  of  ( p )  to  give  M  =  0.2,  0.4,  or  0.6.  The  lower  limit  on  {/>), 
described  earlier,  implies  an  upper  limit  on  the  values  of  M  that  we  can  obtain,  which  for 
our  initial  conditions  is  about  0.6.  Figure  1  shows  the  initial  form  of  the  pressure  field  and 
local  Mach  number,  A4,  for  the  M  =  0.6  case.  It  can  be  seen  that  A4  varies  widely  over 
the  system.  Note  that,  because  of  the  relative  magnitudes  of  the  parameters  in  Eqs.  (6), 
(h)  decreases  as  M  increases  in  our  simulations  according  to 

R  =  _P_  _  2^vl2  (7h\ 

P  ~  ]B|2  7X2!B|2  '  1  } 

By  again  substituting  average  initial  values  for  the  various  parameters,  one  obtains  char¬ 
acteristic  (/?)  values  of  30.0,  7.5,  and  3.3,  respectively. 


C.  Numerical  Algorithm 

In  computations  of  two  dimensional  MHD  systems,  it  is  more  economical  to  replace  the 
magnetic  induction  equation  (Eq.  (lc))  by  the  simpler  magnetic  vector  potential  equation, 

|=vxB  +  ^-V2a,  (lc') 

where  a  is  the  z-component  of  the  magnetic  vector  potential,  and  B  =  V  x  aez.  The 
solenoidality  of  the  magnetic  field  (Eq.  (le))  is  then  trivially  ensured  since  V  •  B  = 
V  •  V  x  aez  =  0.  To  solve  the  governing  equations  (la,  lb,  lc',  and  Id)  our  code, 
CRUNCH2D,  implements  a  Fourier  collocation  method  [42].  The  algorithm  employs  an 
isotropic  truncation  in  Fourier  space  at  each  time-level  [18],  he.,  all  modes  with  |k|  >  y 
are  set  equal  to  zero,  where  k  is  the  Fourier  wavevector  and  N  is  the  number  of  (one  dimen¬ 
sional)  Fourier  modes.  The  modified  Euler  method,  a  second-order  Runge-Kutta  scheme 
[43],  discretizes  time.  The  time  step,  A t,  is  limited  by  a  compressible  MHD  Courant- 
Friedrichs-Lewy  (CFL)  number, 

Ncfl  =  sup{(|Y|+Cf  <  0.3  .  (8) 

where  Cs  =  sound  speed,  and  h  =  2n/N  =  Ax  =  Ay.  Our  simulations  consume  approx¬ 
imately  18  n s  per  time  step  per  grid  point  on  the  NRL  CRAY  X-MP/24.  For  Lundquist 
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numbers  of  50,  100,  and  200,  we  use,  respectively,  642,  1'82,  and  2562  collocation  points. 
The  calculations  with  642  collocation  points  require  approximately  0.12  million  words  of 
core  memory. 

Tire  contrast  between  the  above  spectral  techniques  for  modeling  compressible  flow 
and  finite  difference  algorithms,  which  have  been  in  wider  use  for  compressible  MHD  flows 
[20,  21,  22],  deserves  some  comment.  Finite  difference  algorithms  permit  implementation 
of  complex  boundary  conditions  which  are  presently  more  difficult  or  impossible  for  spec¬ 
tral  methods.  In  addition,  finite  difference  shock-capturing  schemes  are  able  to  maintain 
steeper  and  narrower  shocks  than  their  spectral  counterparts.  Given  that  the  time-scales 
of  interest  in  the  modeling  of  turbulence  are  much  longer  than  shock  transit  times  and 
given  that  we  are  presently  studying  initially  average  subsonic  flows,  this  difference  is  less 
important  than  others. 

On  the  other  hand,  finite  difference  approximations  to  first  and  second  derivatives 
introduce  significant  noise  into  source  terms,  reducing  the  accuracy  of  the  solution.  Solving 
the  magnetic  induction  equation  (Eq.  (lc))  instead  of  that  for  the  vector  potential  (Eq. 
(lc')),  will  ensure  that  only  one  derivative  is  necessary  to  compute  the  current  density, 
at  the  additional  expense  required  to  time-advance  two  dependent  variables  instead  of 
one.  A  careful  formulation  of  the  finite  difference  equations  and  the  differencing  grid 
are  then  necessary  to  maintain  a  solenoidal  magnetic  field  [44].  The  Fourier  collocation 
technique  permits  the  modeling  of  a  much  broader  range  of  spatial  scales  than  do  finite 
difference  techniques  of  comparable  cost.  This  is  a  major  advantage  for  studies  of  both 
compressible  and  incompressible  turbulence  [42],  Finally,  spectral  methods  have  become 
the  method  of  choice  for  general  studies  of  incompressible  turbulence  in  systems  with  simple 
boundary  conditions  [42],  and  the  ability  to  use  similar  algorithms  for  our  compressible 
and  incompressible  calculations  reduces  the  influence  of  numerical  factors  in  comparing 
them. 

A  final  point  regards  the  use  of  the  conservation  form  of  the  model  Eqs.  (1).  The 
mechanical  pressure  is  then  a  “derived”  quantity,  computed  as  the  difference  of  the  total 
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energy  density  and  the  sum  of  the  kinetic  and  magnetic  energy  densities,  as  indicated  by 
Eqs.  (2)  and  (3).  Now  consider  as  examples  (a)  low  /?  regions  in  plasmas,  in  which  the 
mechanical  pressure  is,  by  definition,  much  smaller  than  the  magnetic  pressure,  and  (b) 
local  regions  near  shocks,  where  the  kinetic  energy  can  be  much  larger  than  the  internal 
energy.  In  explicit  schemes,  each  differential  equation  in  the  set,  Eqs.  (1),  is  solved 
separately  from  the  others,  using  parameter  values  from  the  previous  time  step.  If  the 
resulting  solutions  for  the  field  variables  are  not  “synchronized”  in  phase,  errors  in  the 
mechanical  pressure  can  be  as  large  as  the  mechanical  pressure  itself.  We  have  seen  this 
in  finite  difference  calculations,  and  spectral  methods  could  theoretically  have  the  same 
problem.  This  has  driven  some  modelers  to  use  a  pressure  equation  instead  of  Eq.  (Id) 
for  conservation  of  total  energy,  requiring  monitoring  of  the  latter  [45].  A  “fix”  for  this 
problem  in  modeling  shocks  is  the  use  of  sufficient  physical  diffusion  near  the  shocks.  At 
the  Lundquist  and  Mach  numbers  used  here,  the  physical  diffusion  has  been  sufficient  to 
ensure  well-behaved  solutions  for  the  mechanical  pressure. 


D.  Incompressible  Calculations 


To  help  identify  the  turbulent  effects  which  are  due  to  compressibility,  we  run  in¬ 
compressible  calculations  with  a  second  simulation  code.  In  this  limit  the  mass  density 
remains  constant  and  uniform,  and  total  energy  is  not  conserved  in  the  presence  of  finite 
dissipation.  The  form  of  the  incompressible  equations  solved  by  this  second  code  is: 


dv 

dt 


v  x  uiez 


VII  +  jez  x  B  4- 


m 


coupled  with  the  velocity  field  solenoidality  constraint: 


V  v  =  0 


(la') 


and  Eq.  (lc').  The  pressure  head,  II,  is  given  by  p  +  0.5  *  |v|2.  The  integration  of  these 
equations  is  accomplished  with  a  two-dimensional  version  of  the  MHDBOX  explicit  Fourier 
collocation  algorithm  employed  by  Dahlburg  et  al.  [46]  for  studying  the  turbulent  decay  of 


12 


force-free  magnetic  fields.  The  temporal  integration  is  again  modified  Euler.  The  algorithm 
also  employs  an  isotropic  truncation  in  Fourier  space.  We  use  the  same  initial  conditions, 
eqs.  (6),  with  appropriate  adjustments  to  account  for  the  different  normalization. 


III.  RESULTS  AND  DISCUSSION 

Our  investigation  of  compressibility  in  the  Orszag-Tang  vortex  model  relies  primarily 
on  time-dependent  area  integrals,  contour  diagrams,  generalized  streamline  plots,  and 
wavenumber  spectra.  The  temporal  evolution  of  area  integrals  permits  a  direct  comparison 
of  global  compressible  effects  with  their  incompressible  counterparts  and  can  also  show  the 
effects  of  dissipation,  as  parameterized  by  the  Lundquist  numbers.  The  bracket  notation 
represents  a  volume  average,  i.e.,  an  average  of  some  function  (/)  over  the  entire  grid  (Q), 

x,y,t)  dx  dy 

f  dx  dy 
G 

and  thus  applies  to  global  features  of  the  flow.  We  use  contour  plots  and  generalized 
streamline  plots  to  show  how  compressibility  affects  the  local  structures  of  the  flow  field 
and  magnetic  field  and  to  reinforce  inferences  from  the  global  integrals.  The  streamline 
plots  indicate  direction  of  a  vector  function  at  a  given  point,  but  not  the  magnitude  of 
the  function.  Streamline  plots  are  the  best  way  to  visualize  the  form  of  the  compressible 
flow  field.  Since  we  often  compare  the  form  of  the  velocity  field  with  that  of  the  magnetic 
field,  we  also  represent  the  magnetic  field  this  way,  rather  than  using  plots  of  magnetic 
vector  potential  contours.  The  wavenumber  spectra  indicate  interactions  between,  and 
the  importance  of,  fluctuations  at  various  scale  lengths  in  the  magnetofluid.  Since  the 
spreading  out  of  fluctuations  in  Fourier  space  depends  on  the  form  of  the  nonlinearities  in 
the  governing  equations,  we  expect  that  compressibility  will  affect  the  form  of  the  observed 
spectra. 
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A.  Measures  of  Compressible  and  Incompressible  Phenomena 

Some  initial  observations  about  the  global,  integrals  are  necessary.  First  assume  that 
the  appropriate  boundary  conditions  exist  for  surface  integrals  to  vanish.  Because  of  the 
variations  in  mass  density  and  the  dependence  of  the  thermal  pressure  on  temperature 
and  mass  density,  some  quantities  which  are  constants  of  the  motion  in  nondissipative, 
incompressible,  MHD  turbulence  are  no  longer  constants  of  the  motion  in  the  corresponding 
compressible  flows.  In  two  dimensions,  these  quantities  are  (v  •  v  +  B  •  B),  proportional 
to  total  energy  in  an  incompressible  MHD  system  (I?,),  the  mean  square  magnetic  vector 
potential  (a2),  and  the  cross  helicity,  Hc  =  (v  •  B).  Corresponding  conserved  functions  for 
compressible  media  exist  for  the  first  two  quantities:  total  energy,  given  by  Eq.(3),  and  the 
material  integral  (pa2)  [47].  Figure  2  shows  data  on  (pa2)  and  (a2)  for  a  typical  dissipative 
compressible  calculation.  Notice  that  (a2)  fluctuates  about  a  decreasing  monotonic  line 
given  by  (pa2).  The  smooth  decrease  of  the  latter  is  characteristic  of  quantities  that  are 
constant  in  an  ideal  fluid.  The  reason  for  the  proximity  of  the  two  functions  is  the  initial 
condition  p  —  1.  In  general,  the  normalized  ratio  of  inviscid  constants,  (pa2)/(p),  is  best 
for  such  comparisons. 

For  some  boundary  conditions  (e.g.,  doubly  periodic),  Horiuchi  and  Sato  [21]  have 
shown  that  Hc  should  be  an  ideal  constant  of  the  motion  when  the  mass  density  is  uniform 
and  constant  (incompressible  fluid)  or  for  a  homentropic,  polytropic  fluid.  Our  assump¬ 
tion  of  constant  7  is  equivalent  to  the  assumption  that  the  fluid  is  polytropic  (but  not 
homentropic).  In  an  ideal  compressible  fluid,  the  only  sources  of  dissipation  are  shock  dis¬ 
continuities  or  finite  sound  waves,  which  will  steepen  into  shocks  [48].  Since  the  intensity 
of  such  waves  is  smaller  for  lower  Mach  number,  the  fluid  should  have  a  more  spatially 
uniform  entropy  as  Mach  number  decreases,  and  deviations  of  the  pressure  from  the  home¬ 
ntropic,  polytropic  equation  of  state,  Eq.(5),  also  should  be  small.  These  arguments  and 
the  theoretical  arguments  of  Horiuchi  and  Sato,  when  applied  to  the  more  general  class  of 
compressible  media  in  our  calculations,  lead  us  to  expect  the  following  conditions  to  hold: 

(1)  The  plots  of  incompressible  and  compressible  values  of  Hc  versus  time  should  become 
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more  similar  as  Mach  number  decreases  and; 

(2)  The  spatial  variation  of  the  poly  tropic  coefficient,  p/p7,  should  decrease  as  Mach 
number  decreases. 

The  first  statement  indicates  that  Hc  behaves  more  as  an  approximate  ideal  invariant 
as  Mach  number  decreases,  since  we  know  that  it  is  invariant  in  the  incompressible  case. 
Unfortunately  our  results  could  not  rule  out  the  possibility  that  Hc  was  also  invariant  at  the 
higher  Mach  numbers.  The  second  statement  indicates  that,  as  Mach  number  decreases, 
the  fluid  approaches  the  homentropic,  polytropic  case,  for  which  Hc  should  be  invariant. 
The  theory  of  Horiuchi  and  Sato  implies  that,  at  higher  values  of  M,  on  the  other  hand,  Hc 
should  no  longer  be  an  approximate  ideal  invariant  of  the  system.  Further,  the  arguments 
of  Horiuchi  and  Sato,  along  with  our  results,  would  imply  that  the  reason  was  related 
to  the  dissipative  effects  of  compressible  waves  and  the  accompanying  spatial  variation  of 
entropy,  even  in  an  ideal  compressible  fluid  [48].  Both  conditions  (1)  and  (2)  hold  for  our 
calculations.  Figure  3  shows  that  the  evolution  of  cross  helicity  for  a  Mach  number  of 
0.2  is  nearly  identical  to  that  in  the  incompressible  case.  The  plot  for  M  =  0.6  differs 
considerably,  indicating  that  Hc  no  longer  behaves  as  the  ideal,  incompressible  invariant 
does.  From  this  plot,  of  course,  we  cannot  state  with  certainty  that  Hc  is  no  longer  an 
approximate  ideal  invariant  of  the  compressible  system  when  M  is  0.4  or  greater. 

Figure  4  shows  contours  of  the  polytropic  coefficient  for  S  =  100  and  M  =  0.2  and  0.6 
at  early  and  late  times  (<  =  2  and  t  =  8).  This  gives  us  a  picture  of  the  spatially  varying 
entropy  as  well.  For  the  M  —  0.2  case,  the  ratio  of  maximum  and  minimum  values  of 
the  polytropic  coefficient  is  1.27  at  t  =  2,  and  1.12  at  t  =  8.  For  the  M  =  0.6  case, 
the  ratio  of  maximum  and  minimum  values  of  the  polytropic  coefficient  is  about  8.90  at 
t  —  2,  and  2.78  at  t  =  8.  The  ratio  is  more  nearly  one  at  lower  M,  and  thus  the  entropy 
is  more  uniform,  showing  the  approach  to  Eq.  (5)  with  decreasing  Mach  number.  The 
decline  of  this  ratio  with  time  for  both  Mach  number  cases  indicates  that  the  amplitudes 
of  finite  magnetoacoustic  waves  are  decreasing  with  time  as  v  becomes  more  aligned  with 
B  and  as  the  mass  density  becomes  more  uniform.  Apparently  the  system  is  approaching 
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a  quasisteady  state  which  is  also  approximately  homentropic. 

B.  Dynamic  Alignment  and  Selective  Decay 


Dynamic  alignment  is  representative  of  a  more  general  class  of  turbulent  fluid  be¬ 
havior  known  as  self-organization.  Hasegawa  [49]  has  reviewed  the  extensive  work  on  the 
incompressible  MHD  case.  Only  recently  have  studies  of  the  compressible  MHD  case  taken 
place,  using  finite  difference  algorithms  [20,  21,  22]  or  spectral  techniques  [19,  34],  The 
measure  of  correlation  often  used  in  incompressible  calculations  is  the  ratio  of  the  cross 
helicity  ( Hc )  to  incompressible  total  energy  (£’,),  given  by: 

r  He  _  2(v  •  B) 

4  Ei  <|v|2)  +  (|B|2)  • 


This  correlation  coefficient  relates  the  dynamic  alignment  process  to  selective  decay  for 
the  incompressible  case  ( cf .  [50],  [26],  [51]).  In  the  compressible  case,  the  numerator  and 
denominator  which  define  £  are  no  longer  rugged  invariants,  although  analogous  quanti¬ 
ties  may  be  defined.  An  alternative  measure,  which  is  perhaps  more  meaningful  for  the 
compressible  case,  is  the  alignment  factor: 


(vB) 

V(W\2)  <|B|2)  ' 


(11) 


This  is  somewhat  closer  to  the  conventional  definition  of  a  normalized  correlation. 

Figure  5  shows  £  and  a  for  the  incompressible  case  and  as  a  function  of  M  with 
Lundquist  number  S  =  100.  We  see  that,  for  both  the  incompressible  and  compressible 
cases,  v  and  B  become  more  correlated  as  a  function  of  time  and  that  the  growth  of 
correlation  for  M  =  0.2  is  nearly  identical  to  that  of  the  incompressible  system.  As  M 
increases,  however,  the  correlation  at  any  given  time  up  to  t  =  10  is  less.  Some  variation 
from  this  is  seen  in  the  evolution  of  a.  By  both  measures,  though,  compressibility  retards 
the  dynamic  alignment  process  in  the  Orszag-Tang  vortex.  The  variation  with  S  (not 
shown)  is  qualitatively  similar  to  that  found  for  incompressible  calculations  [51].  As  5 
increases,  the  correlation  at  any  given  time  is  less.  Finally,  note  that  a  evolves  in  a  more 
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erratic  manner  than  £,  but  that  both  a  and  £  exhibit  the  same  secular  trend,  viz.,  an 
increase  with  time. 

Past  treatments  of  this  phenomenon  have  attributed  the  temporal  increase  in  align¬ 
ment  in  incompressible  media  to  selective  decay  of  the  energy,  E{  ( ( | v j2 )  +  (|B|2}/2), 
relative  to  the  cross  helicity,  Hc  ((v  •  B))  ([50],  [26],  [51]).  Similar  arguments  apply  to 
the  compressible  case  as  well.  Figure  3  shows,  however,  that  the  functional  dependence 
of  selective  decay  on  M  is  different  from  the  dependence  on  S.  First  consider  variation 
with  5  for  constant  M  and  at  a  particular  time  t.  Figures  3c  and  3d  show  that  both  Hc 
and  E,  increase  with  increasing  5.  However,  the  increase  in  E{  is  greater,  thus  causing  a 
reduction  in  a  and  f  as  5  increases.  This  is  also  the  case  for  incompressible  media  [26]. 

Now  consider  variations  in  Mach  number  in  Figs.  3a  and  3 b.  As  M  increases  with 
constant  S,  Hc  decreases.  The  quantity  Et,  however,  varies  little  with  M.  Figure  5  shows 
that  the  magnetic  and  kinetic  energy  components,  taken  separately,  also  vary  weakly  with 
M.  Because  p  is  0(1),  and  (|v|2)  ~  (|B|2),  the  denominators  of  the  alignment  measures, 
Eqs.  (10)  and  (11),  will  also  show  little  variation  with  Mach  number.  Thus  the  reduction 
in  alignment  with  increasing  M  is  primarily  due  to  a  corresponding  reduction  in  Hc.  The 
discussion  below  of  morphology  cf  the  compressible  flow  reveals  at  least  one  effect  which 
could  cause  this.  As  shown  in  the  next  subsection,  the  velocity  field  exhibits  the  emergence 
of  increased  small  scale  structure  and  greater  structural  variation  with  increasing  M  while 
the  magnetic  field,  in  relative  terms,  does  not.  This  results  in  more  local  reductions  in 
alignment  between  v  and  B  at  a  given  time,  as  M  increases. 

We  also  observe  that  the  dilatation,  V  =  V- v,  is  generally  nonzero  in  the  compressible 
fluid  and  increases  in  magnitude  with  M.  Since  viscous  dissipation  includes  contributions 
from  a  term  dependent  on  the  gradient  of  V ,  the  overall  dissipation  of  the  velocity  field  can 
increase  with  the  amplitude  of  compressible  effects.  This  term  is  important,  for  example, 
in  the  damping  of  acoustic  and  shock  waves.  Such  dissipation  could  lead  to  reductions  in 
(|v)2)  and  Hc  with  increasing  Af.  For  the  total  kinetic  energy  to  remain  approximately 
the  same  as  M  changes  (Fig.  5),  the  change  in  the  local  mass  density  must  compensate 
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for  the  local  changes  in  v. 

C.  Morphology  of  the  Velocity  and  Magnetic  Fields 

The  generalized  streamline  plots  and  the  contour  plots  presented  in  this  section  provide 
a  detailed  visualization  of  the  velocity  and  magnetic  fields  as  the  Orszag-Tang  vortex 
evolves  in  a  compressible  medium.  Again  the  incompressible  solution  provides  a  baseline 
for  interpreting  the  results  for  a  compressible  medium.  Figure  6  shows  the  incompressible 
velocity  field  and  the  corresponding  compressible  velocity  fields  at  t  =  2  as  for  M  = 
0.2,  0.4,  and  0.6  with  S  =  100.  Notice  that  the  M  =  0.2  compressible  solution  closely 
resembles  the  incompressible  solution  in  Fig.  6,  confirming  the  impression  obtained  from 
the  plots  of  global  averages  that  low  Mach  number  MHD  closely  resembles  incompressible 
MHD.  These  two  cases  resemble  each  other  for  all  the  other  times  plotted  (not  shown). 
Also  apparent  in  the  compressible  calculations  is  the  shear  flow  corresponding  to  two  high 
density  jet  structures  emanating  approximately  horizontally  from  the  center  of  the  system, 
at  which  reconnection  is  occurring.  At  M  —  0.6,  the  large  vortex  structures  in  the  middle 
of  the  system  differ  considerably  in  shape  and  smoothness  from  those  present  at  the  lower 
Mach  numbers  and  in  the  incompressible  case.  We  attribute  this  to  differences  in  the  rates 
of  reconnection  occurring  at  the  center  and  edges  of  the  system  for  the  various  calculations 
(next  section).  The  reconnection  and  the  velocity  shear  at  the  center  of  the  system  are 
most  striking  in  Fig.  6  at  the  highest  Mach  number  (M  =  0.6).  Figure  6  also  shows  the 
magnetic  field  for  the  incompressible  and  M  =  0.6  cases  at  the  same  time  ( t  =  2).  This 
field  appears  slightly  less  rounded  as  the  Mach  number  increases,  but  overall  the  magnetic 
field  structure  changes  much  less  than  the  velocity  field  with  increasing  Mach  number. 

Figure  7  compares  the  incompressible  and  M  =  0.2, 0.4,  and  0.6  solutions  for  S  —  100 
at  time  t  =  4,  and  provides  evidence  of  splitting  of  vortices  at  the  higher  Mach  numbers. 
At  lower  Mach  numbers  and  in  the  incompressible  medium,  a  large  vortex  pair  resides 
in  the  middle.  At  M  =  0.6,  each  vortex  has  apparently  split  into  two  separate,  smaller 
vortices  with  an  X-point  in  between.  The  cause  of  this  splitting  with  increased  Mach 
number  apparently  is  the  accompanying  nonlinear  increase  in  the  amplitude  and  steepness 
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of  fluctuations  in  the  fluid  variables,  which  can  cause  appreciable  sources  of  vorticity  at 
smaller  scales  [9].  If  S  is  held  constant  while  M  increases,  the  effects  of  viscous  dissipation 
decrease  relative  to  compressible  sources  of  vorticity. 


We  digress  to  consider  compressible  vorticity  generation  in  greater  detail.  In  a  two- 
dimensional  compressible  magnetofluid,  the  vorticity  equation  is: 


du  _  _  1  Vp  x  Vp  _  __y  1 

—  =  -v  •  Vu>  -  uV  •  v  -f  - - — -  -f  B  •  V-  +  —  V2w 

at  2  p1  p  Sv 


_  .  By.  1  Vp  x  Vp  1 
=  -V  •  (vw  -  — )  +  -  -  +  — V2w 


(12) 


Sv 


In  a  two-dimensional  incompressible  magnetofluid,  the  vorticity  equation  simplifies  to: 


duj  _  _ .  1  _2 

— -  =  -v  •  Vu>  +  B  •  V;  +  — V2w 

C/t  W  V 


(13) 


Unlike  incompressible  MHD,  the  vorticity  quadrupole  production  term  for  compressible 
flow  also  depends  on  the  mass  density  convected  into  the  reconnection  region  (c/.  [52]). 
This  dependence  will  tend  to  enhance  the  quadrupole  over  the  incompressible  result  in 
regions  of  rarefaction  (we  note  parenthetically  that  this  term  By  /p  also  is  significant  for 
dilatation  production,  i.e.,  dV/dt  ex  V  x  By/p).  The  baroclinic  term  will  generate  vor¬ 
ticity  in  any  regions  where  gradients  in  the  mass  density  and  the  pressure  are  not  aligned. 
The  most  prevalent  baroclinic  source  of  vorticity  is  an  asymmetric  pressure  distribution, 
such  as  that  which  arises  at  the  current  sheet  produced  by  reconnection  in  the  center  of 
the  system,  leading  to  larger  amplitude  vorticity  at  the  top  and  bottom  of  the  magnetic 
islands  and  vortex  bifurcation.  Such  pressure  distributions  produce  shocks  or  finite  ampli¬ 
tude  waves,  which  interact  with  each  other  and  with  the  accompanying  nonuniform  mass 
density  distribution  to  produce  or  alter  local  vorticity.  The  high  density  jets  associated 
with  the  reconnection  are  also  a  source  of  significant  velocity  shear  and  vorticity. 

In  addition  to  vortex  splitting,  compressibility  appears  to  generate  more  small  scale 
structure  than  would  otherwise  be  present.  Figure  8  shows  this  by  comparing  the  5  =  200 
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incompressible  and  compressible  M  =  0.6  solutions  at  time  t  =  8.  A  higher  Lundquist 
number,  by  itself,  leads  to  more  structure.  However,  the  M  =  0.6  cases  in  Figures  Id 
and  8b  show  considerable  increases  in  small  scale  structuring  over  their  incompressible 
counterparts.  Figure  9  shows  the  velocity  field  magnitude  for  M  =  0.6  at  t  =  4.  We  see 
that  the  regions  with  the  most  small  scale  vortex  structure  in  the  streamline  plots  also 
have  the  lowest  velocities.  Hence  differences  in  the  small  scale  structure  among  the  various 
cases  do  not  dominate  the  dynamics. 

Figure  10  shows  the  magnetic  field  at  a  later  time,  t  =  8,  as  a  function  of  Mach 
number.  Relative  to  the  velocity  field,  the  magnetic  field  again  shows  the  emergence  of 
much  less  smaller  scale  structure  as  M  increases.  Comparison  with  the  velocity  field  plot 
at  t  =  8  (Fig.  7)  illustrates  the  characteristic  feature  of  dynamic  alignment,  viz.  the 
development  of  a  higher  degree  of  alignment  between  the  magnetic  field  and  the  velocity 
field.  Since  the  velocity  field  shows  more  significant  alteration  in  structure  at  higher  Mach 
numbers,  the  alignment  with  the  magnetic  field  must  be  relatively  less  at  any  given  time. 
This  is  a  possible  explanation  for  the  trend  shown  by  the  global  integrals  in  Fig.  5. 

Figure  11  shows  contours  of  the  dilatation,  P,  for  M  =  0.2  and  0.6  and  5  =  100 
at  t  =  2.  Positive  values  of  P  indicate  regions  of  rarefaction,  negative  values  indicate 
compression.  The  magnitude  of  P  is  a  measure  of  the  level  of  compression  and  rarefaction 
in  the  magnetofluid.  The  dilatation  changes  markedly  as  Af  increases.  The  peak  positive 
and  negative  values,  Pmax  and  |Pmin|,  are  roughly  equal  at  M  =  0.2  while  |Pmin  I  is  greater 
than  Pmax  by  13%  and  55%,  respectively,  at  A1  =  0.2  and  0.4.  In  addition,  Pmax  increases 
by  a  factor  of  5.7  from  M  =  0.2  to  0.4  and  by  a  factor  of  4.3  from  M  =  0.4  to  0.6. 
Similarly  jPm;n  |  increases  by  respective  factors  of  6.2  and  6.0.  These  factors  do  not  give 
an  unambiguous  Mach  number  dependence  for  Pmax  and  Pmin,  but  they  do  show  that  the 
variation  with  M  is  markedly  nonlinear.  |Pmin|  is  usually  greater  than  Pmax  for  M  —  0.6 
at  other  times  t  <  10.  For  the  exceptions,  the  two  are  approximately  equal.  The  trend  of 
increases  of  these  parameters  with  increasing  M  holds  at  all  times  which  we  investigated. 
The  difference  between  |Pmin|  and  Pmax  can  be  as  high  as  a  factor  of  5.  In  addition,  the 
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dilatation  rate  increases  for  increasing  S  and  fixed  M.  For  the  M  =  O.C  case,  as  S  is 
increased  from  100  to  200,  |2?|max  increases  by  a  factor  of  2.5  (at  t  =  2).  Tliis  results 
from  the  accompanying  decrease  in  dissipation  relative  to  the  magnitude  of  compressible 
fluctuations,  thus  reducing  dissipative  smoothing  of  the  velocity  field. 

We  briefly  comment  on  compressible  calculations  with  varying  Lundquist  number  S, 
holding  the  Mach  number  fixed.  The  trends  observed  for  incompressible  calculations  hold 
for  the  compressible  calculations,  as  well.  However,  for  M  =  0.6  and  S  =  200,  the  relative 
amount  of  small  scale  structure  is  higher  than  in  the  corresponding  incompressible  case. 
This  is  evident  in  Fig.  8.  Our  calculations  have  not  covered  sufficiently  large  Lundquist 
numbers  to  determine  whether  the  transition  to  a  fully  developed  turbulent  field  occurs 
earlier  for  higher  values  of  M,  but  the  presence  of  more  small  scale  structure  strongly 
indicates  this  possibility,  possibly  intermittently  in  regions  of  higher  local  Mach  number. 
Note  that  the  plots  of  global  integrals  indicate  that  the  growth  of  correlation  decreases  for 
either  increasing  5  or  M  with  the  other  held  fixed. 

D.  Reconnection  Phenomena 

Figure  12  shows  contour  plots  of  the  electric  current  density  and  the  vorticity  at  time 
t  =  2  for  variable  M  and  with  S  =  100.  The  current  sheet  in  the  center  corresponds  to 
the  sheared  magnetic  field  occurring  there,  as  in  the  incompressible  Orszag-Tang  vortex 
system  [23].  There  is  some  quantitative  variation  in  the  electric  current  density  as  the 
Mach  number  increases,  but  little  qualitative  variation.  This  variation  primarily  reflects  a 
difference  in  the  “scheduling”  of  events  in  the  system.  In  particular,  increased  compress¬ 
ibility  appears  to  retard  the  formation  of  the  current  sheet  centered  on  the  vertical  edges 
of  the  periodic  box.  At  later  times  the  central  current  sheets  stretch  to  the  point  where 
they  exhibit  a  secondary  tearing  instability  [53].  Because  of  the  interaction  between  the 
velocity  and  magnetic  fields  at  this  time,  a  region  of  sheared  flow  coexists  with  the  current 
sheet.  We  see  the  familiar  vortex  quadrupoles  centered  on  the  current  sheet  ([54],  [52]).  In 
addition,  isolated  high  vorticity  regions  form  elsewhere  in  the  box  as  M  increases,  e.g.,  near 
the  upper  left  corner  of  Figure  12d.  The  magnetic  and  kinetic  enstrophics,  ( j 2)  and  (u>2), 
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are  shown  in  Figures  be  and  5/,  respectively,  as  functions  of  time  for  varying  M  and 
5  —  100.  These  averages  show  small  variation  with  Mach  number  and  a  roughly  linear 
variation  with  Lundquist  number  (not  shown  in  the  figure).  As  the  Mach  number  increases 

the  most  notable  features  are  a  decrease  in  the  peak  ( j 2)  and  the  increasing  fluctuations 

/  2  v 

m  ). 

Resistive  heating  in  the  central  region,  due  to  the  current  sheet,  produces  a  massive 
iet  shown  in  Fig.  13  at  t  =  2  for  variable  M,  S  =  100.  As  M  increases:  (1)  the  peak 
Tushy  in  the  jet  increases  (For  the  M  —  0.2  case  peak  density  in  the  jet  is  about  6% 
above  the  initial  value,  while  for  the  M  =  0.6  case  the  peak  density  is  about  110%  above 
the  initial  value);  (2)  the  gradients  in  the  mass  density  increase;  and  (3)  the  jet  becomes 
more  collimated.  It  is  interesting  to  note  that  this  jet  resembles  radio  wavelength  images 
of  astrophvsical  jets,  although  we  do  not  speculate  that  the  cause  of  such  objects  is  the 
.  me.  Figures  13c  and  13d,  at  a  later  time  ( t  =  8),  show  that  the  jet  material  has  spread 
in  the  transverse  direction  and  forms  a  much  less  distinct  structure.  The  temperature  (not 
shown)  contains  evidence  of  both  the  central  vortex  and  magnetic  structures  and  of  the 
jet. 

Also  worthy  of  note  is  the  reconnection  occurring  across  the  left  and  right  (periodic) 
boundaries.  This  reconnection  is  apparent  from  the  magnetic  field  line  plots  in  Figs.  C  and 
10  and  also  from  contour  plots  of  current  density  and  vorticity  in  Fig.  12,  which  show  the 
presence  of  current  and  vortex  sheets.  We  do  not  see  evidence  of  jets  at  the  boundaries, 
:..’u  the  vorticity  plots  do  show  the  existence  of  shear  flows  at  the  “center”  (x  =  7r,  y  =  r) 
and  -edges"  ((0.  j )  and  (2tt,  j /))  of  the  computational  box  As  the  Mach  number  increases  in 
l  :c,  12  and  13,  the  shear  at  the  “edges”  becomes  weaker  relative  to  that  at  the  “center" 
f  the  computational  box. 

F.  Spectral  Distributions 

As  described  above,  compressible  flow  morphology  indicates  the  presence  of  more 
vuail  scale  structure  as  the  initial  average  Mach  number  is  increased.  This  conclusion  also 
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is  borne  out  by  a  study  of  the  volume  averaged  speeiral  distributions  of  various  quantities. 
Two  qxiestions  emerge  naturally:  1)  what  are  the  spectral  distributions  of  compressible 
quantities,  e.g.,  the  mass  density,  since  they  have  no  incompressible  analogue,  and  2)  how- 
does  compressibility  affect  the  spectra  of  quantities  with  incompressible  analogues  [e.g., 
the  magnetic  energy)?  In  general,  at  any  given  time  the  spectral  distributions  show  more 
excitation  in  the  small  spatial  scales  for  all  of  these  quantities  accompanied  by  a  slight 
decrease  in  the  slope  of  the  high- wavenumber  regions  as  M  is  raised.  This  implies  that  the 
usual  caveat  about  increasing  the  numerical  resolution  as  the  Lundquist  (Reynolds-like) 
number  increases  also  should  be  extended  to  the  Mach  number. 

Figure  14  shows  spectral  plots  of  the  mass  density  at  several  different  times  for  the 
M  —  0.2  and  0.6  cases.  One  outstanding  feature  in  these  plots  is  the  rapid  excitation  of  the 
density  at  small  spatial  scales,  implying  that  nonsolenoidal  flows  are  well-developed  after 
one  Alfven  transit  time.  Comparison  shows  that  the  higher  Mach  number  case  for  any 
given  time,  in  a  more  highly  excited  state  (the  level  of  excitation  is  described  roughly  bv 
£p/p  oc  A/2).  Spectra  for  both  Mach  numbers  shown  exhibit  a  rise  in  excitation,  peaking 
at  about  t  =  2,  followed  by  a  gradual  decay. 

The  conclusion  that  increasing  compressibility  leads  to  more  small  scale  structure  is 
straightforward  with  respect  to  a  quantity  such  as  the  mass  density  —  since  the  mass 
density  fluctuations  are  zero  by  hypothesis  in  the  incompressible  case.  What  about  a 
quantity,  like  the  magnetic  field,  which  fluctuates  in  both  cases?  How  will  it  be  affected 
by  variations  in  the  Mach  number?  Figure  15  shows  spectral  plots  of  the  magnetic  energy. 
Here  again  we  find  higher  excitation,  primarily  in  the  short  wavelengths,  as  the  initial 
average  Mach  number  is  increased.  The  longer  wavelengths  are  not  affected  as  strongly,  as 
might  have  been  expected  from  the  field  morphology.  Hence  the  slope  of  the  distribution 
increases  slightly  in  the  high  wavenumber  range  as  the  Mach  number  is  raised,  which 
suggests  that  the  theory  of  Moiseev  et  af.[l6]  might  be  extended  to  the  MHD  case. 

We  next  decompose  the  velocity  field  into  solenoidal  and  nonsolenoidal  parts,  and 
compare  their  spectra  in  Figure  16  at  various  times  for  the  M  =  0.6  case  (c/.[l7]).  To 
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obtain  the  solenoidal  part,  we  compute 

|k  x  v(k)|2 
|k|2  • 

where  the  tilde  denotes  a  Fourier  transformed  quantity.  The  nonsolenoidal  part  is  given 

fcy 

|k-v(k)|2 
|k|2  * 

Our  earlier  conclusion  that  compressible  excitations  develop  quickly  is  borne  out  by  the 
nonsolenoidal  spectra.  In  the  M  —  0.6  case,  there  is  almost  an  equipartition  between 
the  two  velocity  fields,  except  at  the  very  longest  wavelengths.  At  the  later  time  the 
nonsolenoidal  field  is  more  excited  at  the  smaller  scales  than  the  solenoidal  field.  For  the 
M  =  0.2  case  (not  shown),  the  solenoidal  part  of  the  field  generally  is  more  excited  than 
the  nonsolenoidal  part,  especially  in  the  larger  scales. 

A  comparison  with  the  spectral  data  of  Passot  and  Pouquet  (1987)  on  nonmagnetized 
fluids  initialized  with  broadband  noise  can  help  us  to  analyze  our  result  of  an  approximate 
equipartition  between  the  solenoidal  and  nonsolenoidal  components  of  the  velocity  field. 
We  concentrate  on  the  only  case  which  they  presented  for  subsonic  flow,  M  =  0.028, 
although  this  is  far  below  the  range  which  we  investigated.  In  their  simulation,  the  initial 
velocity  field  was  almost  entirely  solenoidal.  For  small  initial  mass  density  fluctuations 
relative  to  M2 ,  Passot  and  Pouquet  found  the  spectrum  to  be  “quasi-incompressible”  at 
all  times.  When  Sp/p  >  M 2  at  t  =  0,  the  initial,  almost  entirely  incompressible,  velocity 
field  changed  in  two  characteristic  times  to  a  state  in  which  the  nonsolenoidal  component 
comprised  over  70%  of  the  kinetic  energy. 

In  our  case,  the  initial  density  fluctuations  are  zero.  In  the  absence  of  the  magnetic 
field,  the  density  would  remain  constant  as  the  kinetic  energy  decayed  to  zero,  and  the 
behavior  should  be  quasi-incompressible,  as  observed  by  Passot  and  Pouquet.  The  dissi¬ 
pation  of  the  magnetic  field,  however,  provides  a  source  of  driving  energy  which  can  cause 
an  appreciable  compressible  component.  The  transient  character  of  this  energy  source 
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and  the  fact  that  the  compressible  system  is  otherwise  identical  to  the  incompressible  sys¬ 
tem  at  very  early  times  and  quite  similar  at  later  times,  indicates  that  an  appreciable 
incompressible  component  could  be  present  as  well.  Under  these  conditions,  approximate 
equipartition  between  the  solenoidal  and  nonsolenoidal  components  seems  reasonable.  At 
the  lower  Mach  numbers,  one  might  expect  a  higher  relative  level  of  solenoidality  than  at 
the  higher  Mach  numbers. 

In  Figure  17  we  show  the  normalized  cross-helicity,  x,  defined  as: 

2  v*(*0- 

iv(*oi2 + iB(*or 

for  M  =  0.2  and  0.6  and  S  =  100  at  early,  middle,  and  late  times.  Incompressible  numerical 
simulations  have  shown  that  x(k)  tends  to  bifurcate  as  time  evolves  ([50],  [26],  [51]).  Since 
\(fc)  oc  Hc(k)/Ej(k ),  this  indicates  that  one  sign  of  cross  helicity  dominates  at  large  spatial 
scales,  while  oppositely  signed  cross  helicity  dominates  at  small  spatial  scales.  We  note 
the  development  of  a  bifurcated  state  (i.e.,  in  which  x  has  different  signs  at  the  large  and 
small  scales)  at  later  times  for  both  Mach  numbers.  For  the  higher  Mach  number  case, 
however,  a  large  central  region  in  the  late  time  plot  shows  no  clear  state  of  correlation. 
This  is  not  found  in  the  lower  Mach  number  plot.  Recall  that  the  level  of  correlation  is  less 
for  the  higher  Mach  number  cases  at  any  given  time  (Figure  5a).  At  t  =  1,  Mach  number 
variation  of  x  is  only  evident  in  the  short  wavelength  regime.  This  regime  is  somewhat 
more  negative  for  the  low  M  run.  However,  at  t  =  5,  the  spectra  exhibit  differences  at  large 
and  small  scales.  The  tendency  toward  bifurcation  is  clear  by  now,  but  domination  by  one 
sign  of  cross  helicity  is  not  as  pronounced  at  short  wavelengths  for  M  =  0.6.  At  t  =  10, 
the  spectra  of  magnetic  and  velocity  fields  in  very  low  and  very  high  wavenumber  regions 
are  correlated,  but  wavenumbers  in  the  middle  range  are  in  a  mixed  state  of  correlation. 
Note  that  total  bifurcation  is  episodic  for  both  cases,  e.g.,  Figure  17a  indicates  that  the 
M  —  0.2  case  has  become  bipolar  at  t  =  5,  however,  at  t  —  10  there  is  no  sharp  break 
between  region  of  oppositely  signed  normalized  cross  helicity.  Numerical  simulations  out 
to  larger  times  are  required  to  see  if  these  systems  settle  into  a  completely  bifurcated  state. 


B  (k) 
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These  simulations  do  show,  however,  that  the  details  of  normalized  cross  helicity  transfer 
are  sensitive  to  the  degree  of  compressibility  of  the  turbulent  magnetofluid. 


IV.  CONCLUSIONS 

To  determine  the  effects  of  compressibility  on  the  phenomenology  of  MHD  turbulence, 
we  have  taken  the  direct  approach  of  choosing  a  simple,  though  nontrivial,  test  problem 
which  has  a  well-known  incompressible  MHD  evolution.  Although  this  choice  of  initial 
conditions  might  limit  the  phase  space  accessible  to  the  compressible  system,  it  does  permit 
a  direct  comparison  of  the  evolution  of  turbulence  in  the  two  cases  and  results  in  the 
identification  of  specific  compressible  effects.  The  present  calculations  have  concentrated 
on  initially  average  subsonic  flows  and  have  permitted  ft  to  change  inversely  with  Mach 
number.  Another  approach  which  we  are  presently  investigating  is  to  maintain  /?  at  a 
constant  value  and  vary  the  Mach  number  by  varying  the  initial  mass  density.  The  use  of 
variable  mass  density  permits  a  wider  range  of  Mach  numbers  to  be  considered,  including 
the  supersonic  regime. 

In  this  paper  we  have  investigated  the  dynamic  evolution  of  compressible  magnetoflu¬ 
ids  as  a  function  of  initial  average  Mach  number,  in  terms  of  global  quantities,  flow  mor¬ 
phology,  and  spectra.  Perhaps  the  most  theoretically  meaningful  quantities  chosen  for 
comparison  are  the  global  quantities.  Comparison  of  global  quantities  shows  that  some 
differences  in  evolution  can  be  described  in  terms  of  initial  average  Mach  number.  These 
results  also  show  that  no  sudden  qualitative  changes  in  behavior  occur  for  the  range  of  ini¬ 
tial  average  Mach  numbers  studied,  but  rather  a  slow  progression.  Our  calculations  thus 
indicate  that,  to  some  extent,  observations  of  compressible  magnetohydrodynamic  tur¬ 
bulence  can  be  interpreted  in  terms  of  incompressible  MHD  turbulence  theory  for  Mach 
numbers  less  than  about  one.  The  relevance  of  incompressible  MHD  turbulence  theory  to 
supersonic  systems  remains  to  be  tested.  Dynamic  alignment  occurs  in  compressible  as 
well  as  in  incompressible  media;  system  evolution  and  the  physics  of  reconnection  are  sim¬ 
ilar,  excluding  variations  in  the  mass  density  and  temperature.  The  alignment  measures 
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of  v  and  B  decrease  at  a  given  time  as  Mach  number  increases,  so  that  compressibility 
retards  or  limits  the  process  for  the  Orszag-Tang  vortex.  The  physics  of  selective  decay 
as  a  function  of  M  is  significantly  different  when  compared  to  the  dependence  on  inverse 
dissipation,  even  though  the  alignment  coefficients  decrease  with  increasing  Mach  number, 
as  with  increasing  Lundquist  number.  Some  variables,  such  as  average  kinetic  and  mag¬ 
netic  energies  and  kinetic  and  magnetic  enstrophies,  show  little  variation  over  the  range 
of  Mach  number  we  have  studied.  Note  that  we  have  drastically  simplified  the  transport 
coefficients  by  assuming  them  to  be  isotropic  and  temperature  independent.  A  more  re¬ 
alistic  trea  tment  of  transport  might  alter  our  conclusions  somewhat.  Note  also  that  the 
behavior  of  ideal,  incompressible  global  invariants  is  not  identical  in  the  two  media.  In  the 
compressible  medium  these  quantities  often  fluctuate  around  the  normalized  compressible 
analogs  or  behave  in  a  manner  that  is  qualitatively  similar  to  the  incompressible  results. 

Our  second  means  of  study  is  point-to-point  local  comparison  of  solutions  as  a  function 
of  Mach  number.  While  the  qualitative  observations  regarding  common  variables  in  the 
Orszag-Tang  vortex,  like  the  velocity  and  magnetic  field,  are  similar,  the  details  of  the 
various  phenomena  differ  in  significant  ways.  Some  features  of  compressible  flows  have  no 
analog  in  incompressible  media.  This  is  true,  for  example,  of  compressible  variations  in 
mass  density.  The  high  density  jets  produced  by  magnetic  reconnection  are  one  important 
feature  which  the  incompressible  calculations  cannot  model  or  estimate.  Another  striking 
feature  in  direct  solution  comparison  is  the  emergence  of  smaller  scale  rotational  structure 
in  the  velocity  field  with  increasing  Mach  number.  The  central  large  scale  vortices  appear  to 
split  into  vortex  pairs  at  high  subsonic  Mach  numbers.  Small  scale  structure  also  increases 
with  higher  Lundquist  numbers  as  in  incompressible  flows  but  more  so  at  higher  M .  The 
size  of  the  Mach  number  relative  to  the  Lundquist  number  should  be  important  in  this 
phenomenon,  since  the  former  is  a  measure  of  the  magnitude  of  velocity  fluctuations  while 
the  latter  is  a  measure  of  the  degree  to  which  dissipation  will  reduce  such  fluctuations.  The 
source  term  in  the  vorticity  evolution  equation  for  compressible  flow  (Eq.  (12))  contributes 
to  changes  in  the  vorticity  through  pressure  and  density  fluctuations  [9],  as  well  as  by 
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altering  the  magnetic  source  term.  These  changes  should  translate  into  smaller  scale 
structure  in  the  flow  field.  Spectral  plots  as  a  function  of  M  support  this  conclusion  and 
hint  that  compressibility  might  significantly  alter  spectral  transfer  in  the  inertial  range. 

None  of  the  above  results  have  depended  explicitly  upon  the  feature  usually  associ¬ 
ated  with  compressible  flows — finite  sound  waves  and  shock  waves.  We  saw  little  direct 
evidence  of  such  waves  in  our  calculations.  This  will  be  the  subject  of  our  next  paper  on 
the  supersonic  regime.  Our  earlier  results  [19],  combined  with  those  of  Passot  and  Pou- 
quet,  indicate  that  novel  phenomena  occur  when  compressible  initial  conditions  are  used. 
However,  we  point  out  that,  in  a  region  containing  shocks  and  for  which  (v)  =  0,  the  flows 
will  be  subsonic  in  the  portions  bounded  by  shocks,  so  that  the  conclusions  that  we  have 
reached  here  will  be  important  to  the  case  of  supersonic  flow  as  well. 
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Figure  1.  Initial  conditions:  a.  velocity  field;  b.  magnetic  field;  Initial  conditions  specific 
to  the  M  =  0.6  case:  c.  thermal  pressure  field;  d.  local  Mach  number.  Minimum 
and  maximum  values  are  0.73  and  7.9  for  thermal  pressure  and  0.0  and  0.97  for  local 


Mach  number. 


iP 

tX) 


Figure  2.  Time  evolution  of  the  mean  square  magnetic  vector  potential  ((a2)  -  solid 
line)  and  its  compressible  analogue  ({pa?)  -  dashed  line),  both  for  the  compressible 
M  =  0.6,  S  —  100  run. 
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Figure  3.  Variation  of  time  evolution  of  global  quantities  parameterized  by  Mach  number 
(M)  and  by  Lundquist  number  (5):  a.  variation  of  cross  helicity  with  M  (5  =  100) 
(solid  line— incompressible  result,  dashed  line- M  =  .2,  mixed  dashed  Yme-M  —  .4,  and 
dash-dotted  line-M  =  .6);  b.  variation  of  fluctuating  energy  (magnetic  +  kinetic)  with 
M;  (figure  key  same  as  a);  c.  variation  of  cross  helicity  with  5  ( M  =  .6)  (triangles- 
S  =  200,  circles-S  =  100,  and  squares-5  =  50);  and  d.  variation  of  fluctuating  energy 
(magnetic  +  kinetic)  with  S  (figure  key  same  as  c). 
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Figure  5.  Variation  of  time  evolution  of  global  quantities  parameterized  by  Mach  number 
(A/),  S  =  100  for  all  cases:  (solid  line-incompressible  result,  dashed  line-A/  =  .2, 
mixed  dashed  line—  M  =  .4,  and  dash-dotted  line- -M  =  .6);  a.  correlation  coefficient; 
b.  alignment  factor;  c.  kinetic  energy;  d.  magnetic  energy;  e.  kinetic  enstrophy;  and 
f.  magnetic  enstrophy. 
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Figure  6.  Variation  of  fields  with  Mach  number  (M)  at  t  =  2  (5  =  100):  a.  velocity 

field  (incompressible  case);  b.  velocity  field  (M  =  0.2);  c.  velocity  field  (M  =  0.4); 
d.  velocity  field  (M  =  0.6);  e.  magnetic  field  (incompressible  case);  and  f.  magnetic 
field  {M  =  0.6). 


Figure  7.  Variation  of  velocity  field  with  Mach  number  ( M)  at  t 
incompressible  case;  b.  M  =  0.2;  c.  M  -  0.4;  and  d.  M  —  0.0. 
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c.  d. 

Figure  12.  Variation  of  fields  at  t  =  2  (5  =  100)  a.  electric  current  density  (incompressible 
case);  b.  electric  current  density  (M  =  0.6);  c.  vorticity  (incompressible  case);  and  d. 
vorticity  ( M  =  0.6).  Minimum  and  maximum  values  are  a.  approximately  -20.0  and 
20.0,  b.  -13.1  and  20.4,  c.  approximately  -6.40  and  6.40,  and  d.  -8.10  and  6.45. 
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Figure  13.  Variation  of  the  mass  density  as  a  function  of  Mach  number  M  ( 5  —  100) 
t  =  2  (M  =  0.2);  b.  t  =  2  (M  =  0.6);  c.  t  -  8  (M  =  0.2);  and  d.  <  =  8  (M  =  0.6). 


Figure  16.  Spectral  plots  of  the  solenoidal  and  nonsolenoidal  kinetic  energy  density  for  the 
M  —  0.6  case:  solenoidal  part  at  a.  t  =  2;  b.  t  =  4;  c.  t  =  6;  d.  t  =  8;  Nonsolenoidal 
part  at  e.  t  =  2;  f.  t  =  4;  g.  t  =  6;  h.  t  =  8.  S  =  100  for  all  cases.  The  dashed  line 
has  a  reference  slope  of  -4/3. 
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Figure  17.  Spectral  plots  of  the  normalized  cross-helicity:  M  =  0.2  case  at  a.  t  =  1;  b. 
t  —  5;  c.  t  =  10;  M  =  0.6  case  at  d.  t  —  1;  e.  t  =  5;  f.  i  —  10.  S  =  100  for  all  cases. 
The  dashed  line  has  the  reference  value  zero. 
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